library(Seurat)
## Attaching SeuratObject
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.0 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.1.8
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(magrittr)
##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extract
library(liana)
library(cowplot)
##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamp
library(BiocParallel)
library(nichenetr)
Load Rds with NK cells, liana results and filter the interactions by ICR
dat <- readRDS("Data/MM_atlas.Rds")# scRNA-Seq atlas
receiver = c("NK exhausted", "NK resident")
## labels
library(readr)
all_labels_df <- read_csv("all_labels_df.csv",
col_types = cols(...1 = col_skip()),
trim_ws = FALSE)
## New names:
## • `` -> `...1`
all_labels_df <- as.data.frame(all_labels_df)
rownames(all_labels_df) <- all_labels_df[,"keyName"]
all_labels_df$keyName <- NULL
colnames(all_labels_df) <- "labels_interactions"
all_labels_df$labels_interactions <- as.factor(all_labels_df$labels_interactions)
dat <- AddMetaData(dat, all_labels_df)
dat <- SetIdent(dat, value = "labels_interactions")
ICR = c("ADCY7|ADGRE5|ADGRG1|BTN3A2|CD244|CD47|CD53|CD81|CD96|CLEC2D|CXCR4|DIP2A|FAS|HAVCR2|HLA-DPB1|IGF2R|KIR2DL1|KIR2DL3|KLRB1|KLRC1|KLRC1_KLRD1|KLRF1|KLRG1|LAG3|LAIR1|LILRB1|PTPRC|RACK1|SIGIRR|TIGIT|TNFRSF14")
interactions <- read_csv("all_interactions_defaults.csv", col_types = cols(...1 = col_skip())) %>%
filter(target %in% receiver) %>%
filter(grepl(ICR,receptor.complex)) %>%
filter(!grepl("KLRC2|KLRC3", receptor.complex))
## New names:
## • `` -> `...1`
Load models for NicheNet:
ligand_target_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds"))
ligand_tf_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_tf_matrix_nsga2r_final.rds"))
lr_network = readRDS(url("https://zenodo.org/record/7074291/files/lr_network_human_21122021.rds"))
sig_network = readRDS(url("https://zenodo.org/record/7074291/files/signaling_network_human_21122021.rds"))
gr_network = readRDS(url("https://zenodo.org/record/7074291/files/gr_network_human_21122021.rds"))
weighted_networks = readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final.rds"))
weighted_networks_lr = weighted_networks$lr_sig %>% inner_join(lr_network %>% distinct(from,to), by = c("from","to"))
# Select ligands actively targetting ICR
ligands = unique(sort(interactions$ligand.complex))
ligands = ligands[ligands %in% colnames(ligand_target_matrix)]
ligands
## [1] "B2M" "CALM1" "CCN2" "CD22" "CD40LG" "CD48"
## [7] "CD55" "CD99" "CDH1" "CEACAM1" "CLEC2B" "CLEC2D"
## [13] "CXCL12" "FAM3C" "GNAS" "GZMB" "HLA-A" "HLA-B"
## [19] "HLA-C" "HLA-DPB1" "HLA-DQA1" "HLA-DQA2" "HLA-DQB1" "HLA-DRA"
## [25] "HLA-DRB1" "HLA-DRB5" "HLA-E" "HMGB1" "IL1B" "LGALS1"
## [31] "LGALS3" "LGALS9" "LILRB4" "LTF" "MIF" "NECTIN1"
## [37] "NECTIN2" "NRG1" "SIRPA" "TGFB1" "TGM2" "THBS1"
## [43] "TNFSF13" "TNFSF13B"
dat_nk = subset(dat, subset = labels_interactions %in% receiver)
I will consider a gene expressed if it is present in at least 10% of cells for each celltype as suggested in NicheNet Documentation
## receiver
expressed_genes_receiver = get_expressed_genes(receiver, dat_nk, pct = 0.10)
background_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
As target genes I am going to use the genes included in those regulons that contains ICRs mediators from the TF analysis. This is aimed to integrate the two endpoints, signalling and regulatory network, to buil the signalling cascade of ICRs:
geneset_oi <- read_tsv("Data/all_ICR-mediators_target_genes.csv", col_types = cols(), col_names = "gene") %>%
pull(gene) %>%
.[. %in% rownames(ligand_target_matrix)]
nichenet_activities <- predict_ligand_activities(
geneset = geneset_oi,
background_expressed_genes = background_genes,
ligand_target_matrix = ligand_target_matrix, potential_ligands = ligands
)
# prepare data for visualization
vis_liana_nichenet <- interactions %>%
inner_join(nichenet_activities, by = c("ligand.complex" = "test_ligand")) %>%
arrange(pearson) %>%
mutate(ligand = fct_inorder(ligand.complex))
# prepare NicheNet figure
nichenet_scores_plot <- vis_liana_nichenet %>%
group_by(ligand) %>%
summarize(pearson = mean(pearson)) %>%
ggplot(aes(y = ligand, x = pearson)) +
geom_bar(stat = "identity") +
ggtitle("NicheNet") +
xlab("Pearson's score") +
theme_cowplot() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
axis.line.y = element_line(color = "white"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))
# prepare LIANA figure
liana_receptor_heatmap <- vis_liana_nichenet %>%
ggplot(aes(y = ligand, x = receptor.complex, fill = aggregate_rank)) +
geom_tile() +
theme_cowplot() +
ggtitle("LIANA") +
ylab("Ligand") + xlab("Receptor") +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1),
plot.title = element_text(hjust = 0.5),
panel.grid.major = element_line(colour = "gray", linetype = 2),
legend.position = "left")
# combine plots
plot_grid(liana_receptor_heatmap, nichenet_scores_plot,
align = "h", nrow = 1, rel_widths = c(0.8,0.3))
nichenet_activities %>% arrange(-pearson)
## # A tibble: 44 × 5
## test_ligand auroc aupr aupr_corrected pearson
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 THBS1 0.657 0.0634 0.0566 0.134
## 2 LGALS3 0.589 0.0546 0.0478 0.0976
## 3 CD40LG 0.707 0.0309 0.0241 0.0876
## 4 CCN2 0.630 0.0530 0.0462 0.0717
## 5 TNFSF13 0.629 0.0201 0.0133 0.0708
## 6 HLA-DRB5 0.634 0.0260 0.0192 0.0666
## 7 CD22 0.666 0.0191 0.0123 0.0629
## 8 HLA-C 0.681 0.0175 0.0107 0.0601
## 9 LGALS9 0.666 0.0164 0.00958 0.0552
## 10 SIRPA 0.667 0.0260 0.0192 0.0544
## # … with 34 more rows
best_upstream_ligands = nichenet_activities %>% arrange(-pearson) %>% pull(test_ligand)
active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()
active_ligand_target_links_df = na.omit(active_ligand_target_links_df)
active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.25)
nrow(active_ligand_target_links_df)
## [1] 100
head(active_ligand_target_links_df)
## # A tibble: 6 × 3
## ligand target weight
## <chr> <chr> <dbl>
## 1 THBS1 BCL2L1 0.0152
## 2 THBS1 PECAM1 0.0668
## 3 THBS1 SOCS3 0.0117
## 4 LGALS3 BCL2L1 0.0102
## 5 LGALS3 SOCS3 0.0649
## 6 CD40LG BCL2L1 0.0252
order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev()
order_targets = active_ligand_target_links_df$target %>% unique()
vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()
p_ligand_target_network = vis_ligand_target %>% make_heatmap_ggplot("Microenvironment ligands","Inhibitory target genes", color = "purple",legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") + scale_fill_gradient2(low = "whitesmoke", high = "purple", breaks = c(0,0.005,0.01)) + theme(axis.text.x = element_text(face = "italic"))
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
p_ligand_target_network
# get the ligand-receptor network of the top-ranked ligands
receptors = lr_network %>% pull(to) %>% unique()
expressed_receptors = intersect(receptors,background_genes)
lr_network_top = lr_network %>% filter(from %in% best_upstream_ligands & to %in% expressed_receptors) %>% distinct(from,to)
best_upstream_receptors = lr_network_top %>% pull(to) %>% unique()
# get the weights of the ligand-receptor interactions as used in the NicheNet model
lr_network_top_df = weighted_networks$lr_sig %>% filter(from %in% best_upstream_ligands & to %in% best_upstream_receptors)
# convert to a matrix
lr_network_top_df = lr_network_top_df %>% spread("from","weight",fill = 0)
lr_network_top_matrix = lr_network_top_df %>% select(-to) %>% as.matrix() %>% magrittr::set_rownames(lr_network_top_df$to)
# perform hierarchical clustering to order the ligands and receptors
dist_receptors = dist(lr_network_top_matrix, method = "binary")
hclust_receptors = hclust(dist_receptors, method = "ward.D2")
order_receptors = hclust_receptors$labels[hclust_receptors$order]
dist_ligands = dist(lr_network_top_matrix %>% t(), method = "binary")
hclust_ligands = hclust(dist_ligands, method = "ward.D2")
order_ligands_receptor = hclust_ligands$labels[hclust_ligands$order]
vis_ligand_receptor_network = lr_network_top_matrix[order_receptors, order_ligands_receptor]
p_ligand_receptor_network = vis_ligand_receptor_network %>% t() %>% make_heatmap_ggplot("Microenvironment ligands","Inhibitory receptors", color = "mediumvioletred", x_axis_position = "top",legend_title = "Prior interaction potential")
p_ligand_receptor_network
active_signaling_network = get_ligand_signaling_path(ligand_tf_matrix = ligand_tf_matrix, ligands_all = ligands, targets_all = geneset_oi, weighted_networks = weighted_networks)
data_source_network = infer_supporting_datasources(signaling_graph_list = active_signaling_network,lr_network = lr_network, sig_network = sig_network, gr_network = gr_network)
head(data_source_network)
## # A tibble: 6 × 5
## from to source database layer
## <chr> <chr> <chr> <chr> <chr>
## 1 ABCA1 CRK harmonizome_GEO_GENE harmonizome_gr regulatory
## 2 ABCA1 DGKZ harmonizome_GEO_GENE harmonizome_gr regulatory
## 3 ABCA1 INPP5D harmonizome_GEO_GENE harmonizome_gr regulatory
## 4 ACHE CRK harmonizome_GEO_GENE harmonizome_gr regulatory
## 5 ACHE PIK3CA harmonizome_GEO_GENE harmonizome_gr regulatory
## 6 ACHE PTPN11 harmonizome_GEO_GENE harmonizome_gr regulatory
output_path = "Results/all_NK_cells/"
dir.create(output_path, recursive = TRUE)
## Warning in dir.create(output_path, recursive = TRUE): 'Results/all_NK_cells'
## already exists
write_output = TRUE # change to TRUE for writing output
# weighted networks ('import network' in Cytoscape)
save_nets <- function(active_signaling_network, ligands, targets, output_path, write_output){
if(write_output){
bind_rows(active_signaling_network$sig %>% mutate(layer = "signaling"), active_signaling_network$gr %>% mutate(layer = "regulatory")) %>% write_tsv(paste0(output_path,"weighted_signaling_network.txt"))
}
# networks with information of supporting data sources ('import network' in Cytoscape)
if(write_output){
data_source_network %>% write_tsv(paste0(output_path,"data_source_network.txt"))
}
# Node annotation table ('import table' in Cytoscape)
specific_annotation_tbl = bind_rows(
tibble(gene = ligands, annotation = "ligand"),
tibble(gene = targets, annotation = "target"),
tibble(gene = c(data_source_network$from, data_source_network$to) %>% unique() %>% setdiff(c(targets,ligands)) %>% intersect(lr_network$to %>% unique()), annotation = "receptor"),
tibble(gene = c(data_source_network$from, data_source_network$to) %>% unique() %>% setdiff(c(targets,ligands)) %>% intersect(gr_network$from %>% unique()) %>% setdiff(c(data_source_network$from, data_source_network$to) %>% unique() %>% intersect(lr_network$to %>% unique())),annotation = "transcriptional regulator")
)
non_specific_annotation_tbl = tibble(gene = c(data_source_network$from, data_source_network$to) %>% unique() %>% setdiff(specific_annotation_tbl$gene), annotation = "signaling mediator")
if(write_output){
bind_rows(specific_annotation_tbl,non_specific_annotation_tbl) %>% write_tsv(paste0(output_path,"annotation_table.txt"))
}
}
sub_nets <- function(ligand_tf_matrix,ligand, targets){
active_signaling_network = get_ligand_signaling_path(ligand_tf_matrix = ligand_tf_matrix,
ligands_all = ligand,
targets_all = targets,
weighted_networks = weighted_networks)
data_source_network = infer_supporting_datasources(signaling_graph_list = active_signaling_network,
lr_network = lr_network,
sig_network = sig_network,
gr_network = gr_network)
active_signaling_network_min_max = active_signaling_network
active_signaling_network_min_max$sig = active_signaling_network_min_max$sig %>% mutate(weight = ((weight-min(weight))/(max(weight)-min(weight))) + 0.75)
active_signaling_network_min_max$gr = active_signaling_network_min_max$gr %>% mutate(weight = ((weight-min(weight))/(max(weight)-min(weight))) + 0.75)
graph_min_max = diagrammer_format_signaling_graph(signaling_graph_list = active_signaling_network_min_max, ligands_all = ligand, targets_all = targets, sig_color = "indianred", gr_color = "steelblue")
# To render the graph: uncomment following line of code
graph = DiagrammeR::render_graph(graph_min_max, layout = "kk", title = ligand)
output_path = paste0("Results/all_NK_cells/", ligand, "/")
dir.create(output_path, recursive = TRUE)
write_output = TRUE # change to TRUE for writing output
save_nets(active_signaling_network,ligands = ligand, targets = targets, output_path, write_output)
return(graph)
}
plot_list <- htmltools::tagList()
for (l in 1:length(ligands)){
graph = sub_nets(ligand_tf_matrix, ligands[l], geneset_oi)
plot_list[[l]] <- graph
}
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## ℹ The deprecated feature was likely used in the DiagrammeR package.
## Please report the issue at
## <]8;;https://github.com/rich-iannone/DiagrammeR/issueshttps://github.com/rich-iannone/DiagrammeR/issues]8;;>.
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/B2M' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CALM1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CCN2' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD22' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD40LG' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD48' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD55' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD99' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CDH1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CEACAM1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CLEC2B' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CLEC2D' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CXCL12' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/FAM3C' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/GNAS' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/GZMB' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-A' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-B' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-C' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DPB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DQA1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DQA2' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DQB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DRA' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DRB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DRB5' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-E' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HMGB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/IL1B' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LGALS1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LGALS3' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LGALS9' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LILRB4' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LTF' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/MIF' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/NECTIN1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/NECTIN2' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/NRG1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/SIRPA' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/TGFB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/TGM2' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/THBS1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/TNFSF13' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/TNFSF13B' already exists
plot_list
The interactions with LAG3, TIGIT, LAIR1 and other ICR looks interesting. Because of problems visualizing the whole network on R I’m going to continue the analysis on python.
I built the sub-networks with all the target genes (all the genes included in those regulons that were including ICRs mediators). This was more for testing and visualization purposes. Selecting specific ligands or, idndirectly receptors, and some target genes (e.g. SHP1/SHP2), maybe we could identify smaller sub-networks with clearer pathways. Now im going to focus on the network dedicated to only exh NK cells in MM and check which ligands and target genes could better suit our aims.